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I.  INTRODUCTION 


Missile  simulation  is  an  indispensable  tool  with  which  the  engineer’s 
computer  can  emulate  the  behavior  of  an  actual  missile  in  flight.  Because  a 
simulation  is  usually  not  an  exact  reproduction  of  a  physical  system,  the  engi 
neer  must  make  certain  assumptions  or  approximations  as  to  how  a  physical  sys¬ 
tem  should  be  mathematically  represented  in  a  computer  program.  The  results 
of  the  simulation  should  allow  the  engineer  to  verify  the  performance  of  the 
system  and  to  validate  his/her  assumptions. 

The  most  straightforward  method  of  tackling  a  missile  system  simulation 
would  be  to  divide  the  system  into  discrete  sections  and  construct  subrou¬ 
tines  that  would  accurately  model  their  tangible  equivalent.  The  subroutines 
that  wil]  be  addressed  in  this  report  are  the  fin  actuators  and  the  numerical 
integration  scheme. 

II.  FIN  ACTUATOR  MODELS 

Two  types  of  fin  actuator  models  that  we  will  be  discussing  are  a  linear 
second-order  model  and  a  non-linear  second-order  model. 


A.  Linear  Model 


A  simple  model  that  could  describe  an  actuator's  dynamics  is  a 
linear  second-order  system  with  damping  zeta  (5)  and  natural  frequency  omega 
( 0*1)  -  Such  a  system  exemplifies  unity  gain  out  to  ujq,  peaking  of  gain  at 
inversely  proportional  to  damping,  a  -20dB  per  decade  slope  in  gain  after 
and  a  180°  change  in  phase.  Step  response  characteristics  such  as  rise 
time  and  overshoot  are  a  function  of  bandwidth  (ton)  and  damping. 

The  transfer  function  of  a  second-order  system  is  given  below,  where 
6  is  the  output  and  <$<.  is  the  input.  Figure  1  shows  one  of  many  possible 
methods  of  implementing  the  transfer  function  as  a  block  diagram. 
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Figure  1.  Second-order  block  diagram. 
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The  differential  equation  describing  the  system’s  dynamics  is: 


5  -  u*i2  (  &c  ~  6  ~  6 — 

m 

The  differential  equation  when  encoded  in  FORTRAN  becomes: 

DEFDD  *  (OMEGA**2)*(DEFC-D£F-(DEFD*2*ZETA/OMEGA) ) 

The  system's  state  equations  are  as  follows: 
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Where  6  (xi)  is  the  fin  position,  6  (x2)  is  the  fin  velocity,  and 
6  (x2)  is  the  fin  acceleration. 

1.  Analytical  Solution 


Since  there  will  be  differences  between  the  solutions  to  the  dif¬ 
ferential  equations  for  the  actual  and  simulated  systems,  we  will  develop  an 
analytical  solution  for  comparison.  The  input  (6C)  will  be  a  unit  step. 


5c(c)  *  1 

* 

«c(S) 


from  the  transfer  function: 


6(S)  -  G(S)*5c(S) 


S(S2  +  2 £u>n  S  +  oft2) 


6(t)  -  L~i[  6(S)  ]  -. 


5(t): 


From  a  table  of  Laplace  transform  pairs,  we  find  the  solution  to 


<5(  t ) 


<%2  <%2  Vi  -  c2 


e  sin(  ain  C2  +  cos  ^ 
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1  m  . 

6(t)  -  1  -  1  “  c2  +  cos  5) 

Vi-  c2 
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For  a  solution  to  the  first  integral,  6(t): 


assume  6(0+)  =  0. 


6(0  -  L_1[S5(S)  -  6(0+)] 


S  6(S) 


<4n 


+  2  5o)nS  +  ton^ 


get: 


From  the  Laplace  transform  pairs  table: 


*6(t )  =  ufe2 


*6(0 


^  V1- c2 

mie  ca*lt  sin((JftVl-g2  t) 

Vl-'T2 


e"^C  si n(ain  V1"?2  O 


Using  an  of  144  rads/sec,  and  a  i;  of  0.6,  as  an  example  we 


6(t)  =*  1-1.25  e~86*4t  sin(  1 15 . 2t  +  0.9273)  rad 
6(t)  *  180  e  88'4t  sin(115.2t)  rad/sec 
2.  Simulation  Results 


Output  data  f’-om  the  simulation  were  tabulated  and  plotted  to 
produce  the  graphs  shown  in  Figures  2,  3,  and  4.  Figure  2  shows  fin  position 
(angle)  versus  time.  Figure  3  is  a  plot  of  fin  velocity  versus  time  and  Fig¬ 
ure  4  shows  fin  acceleration  versus  time.  These  outputs  are  in  response  to  a 
step  input.  Comparison  of  analytical  and  simulated  solutions  are  given  in 
Section  III.  The  FORTRAN  code  for  this  program  and  a  short  discussion  are 
given  in  the  Appendix. 

B.  Non-Linear  Model 


A  second  type  of  model  is  one  that  contains  physical  limitations, 
which  were  added  to  the  linear  second-order  model  to  yield  a  second-order  non¬ 
linear  model. 


1.  Development 

The  second-order  linear  model  was  modified  to  include  character¬ 
istics  typical  of  an  actuator  motor.  This  resultant  non-linear  model  more 
closely  emulates  the  real  thing.  These  characteristics  are  inherent  limita¬ 
tions  of  the  physical  system  and  are  non-linear.  They  include;  position 
limits  (fin  stops),  velocity  limits  (slew  rate  limits),  acceleration  limits 
(finite  torque),  a  deadband  in  the  rate  feedback  (hysteresis),  and  aerodynamic 
hinge  moments.  Figure  5  shows  a  block  diagram  with  the  note  describing  the 
differential  equation. 
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of  fin  angle  versus  time  for  the  linear  mode  1 . 
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Figure  3.  Plot  of  fin  velocity  versus  time,  linear  model. 
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Figure  5.  Second  order  non-linear  block  diagram. 

NOTE:  Descriptive  differential  equation  is: 

5  *  Oh2  ( 6c  -  Slim  -  Ratefb  -  6  -  mi 

The  describing  differential  equation  that  appears  in  the  simula¬ 
tion  code  is: 

DEFDD  =  ( 0WEGA**2 ) * (DEFC  -  DEF  -  RATEFB  -  (DEFD*2 . *ZETA/ OMEGA) )  -  HM 
where:  DEF  =*  <5  =  fin  position 

DEFD  =  6  =  fin  velocity 
DEFDD  =*  <5  =»  fin  acceleration 
HM  *  hinge  moment 
RATEFB  *  rate  feedback 
OMEGA  -  bandwidth  (  ^ 

ZETA  *  damping  factor. 

As  can  be  seen  from  the  block,  diagram,  no  trivial  analytical 
solution  for  the  differential  equation  exists,  so  we  must  rely  on  a  numeri¬ 
cally  integrated  solution. 
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2.  Simulation  results 


The  fin  angle,  velocity,  and  acceleration  versus  time  for  the 
non-linear  second-order  model  were  plotted  and  appear  in  Figures  6,  7,  and  8. 
The  hinge  moment  for  these  plots  was  set  to  zero. 

Compared  to  the  linear  model,  the  position  rise  time  is 
lengthened  due  to  the  rate  and  acceleration  limits  of  5.25  rads/sec  and  300 
rads/sec^,  respectively.  The  overshoot  is  less  due  to  the  position  limit  of 
0.436  rads.  Plots  for  the  negative  fin  deflection  were  not  given  since  the 
results  were  a  mirror  image  of  the  positive  deflection. 

Various  hinge  moment  constants  were  added  and  the  results  are 
shown  in  Figures  9  through  12.  The  hinge  moment  in  Figure  9  is  an  aiding 
moment:  Note  the  increased  overshoot  and  faster  rise  time.  The  hinge  moments 

in  Figures  10  through  12  are  hindering  moments:  Note  that  some  of  the  plots 
show  that  the  commanded  fin  deflections  were  not  achieved  and  in  Figure  12  a 
limit  cycle  resulted.  The  complete  FORTRAN  code  for  this  non-linear  model  is 
included  in  the  Appendix. 
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Figure  6.  Plot  of  fin  position  vs.  time,  non-linear. 


RADS/SEC  (*E  +01) 


Figure  7.  Plo 


11 


ANGLE  (*E  -00) 


CD 

CO 

S3 


CD 

in 

S) 

LH 

S3 

m 

co 

s 

S3 

S3 

S3 

S3 

S3 

S3 

12 


Figure  9.  Plot  of  fin  angle  vs.  time,  hinge  moment  k  =  -800. 
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Figure  10.  Plot  of  fin  angle  vs.  time,  hinge  moment 
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III.  NUMERICAL  INTEGRATION 


A  method  of  solving  differential  equations  in  a  computer  program  is  to 
find  the  solutions  by  numerical  integration.  This  process  usually  involves 
iterative  mathematical  calculations. 

A.  Development 

Solutions  to  mathematically  modeled  systems  with  transient  or  time- 
dependant  processes  usually  involves  solving  differential  equations.  Some 
differential  equations  can  be  solved  analytically,  however  most  cannot,  par¬ 
ticularly  if  non-linear  elements  are  present.  Therefore,  in  computer 
programs,  numerical  Integration  methods  can  be  used  to  solve  differential 
equations  with  acceptable  accuracy. 

One  numerical  integration  method  commonly  used  is  the  Runge-Kutta 
fourth-order  method  (R-K  4).  There  are  many  good  texts  available  that  inves¬ 
tigate  and  develop  this  method  in  detail,  however,  only  a  general  overview 
will  be  given  here. 

The  R-K  4  scheme  uses  a  weighted  average  of  four  estimates  to  calcu 
late  an  approximation  to  the  solution.  A  simple  form  of  a  first  order  differ¬ 
ential  equation  is: 


dx 

dt 


f(x,t) 


» 


where  f(x,t)  is  a  known  function.  Substituting  values  xn  and  tn  into  the 
equation  we  get  the  slope  of  the  solution  curve  at  a  known  starting  point,  tn 

The  R-K  4  approach  is  to  find  an  approximate  slope  at  a  known 
starting  (xn,tn)  and  use  this  slope  and  a  small  time  increment  to  proceed  to 
the  next  point.  Then  assuming  that  this  new  point  is  a  known  point  on  the 
solution  line,  again  find  an  approximate  slope  at  this  new  point  and  proceed 
to  the  next  point  by  incrementing  the  time  step.  The  equations  are: 

k^  *  f(tn,xn) 

h  klh 

k£  *  f(tn  +  j,  xn  +  — j-) 

h  k2h 

k3  -  f(tn  +  y,  xn  +  -y— ) 

k4  *  f(tn  +  h,  xn  +  k3h) 

xn+1  *  xn  +  -jr  (k^  +  2k£  +  2k3  +  k4) 


where:  h  is  the  time  increment 
tn  is  the  initial  time 
xn  is  the  initial  solution  point. 
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B.  Accuracy 


The  accuracy  of  the  solution  is  a  function  of  the  time  step,  h.  A 
smaller  time  step  results  in  better  accuracy  but  the  program  run  times  are 
increased  compared  to  a  larger  time  step.  Too  small  a  time  step  can  cause 
round-off  or  truncation  errors  due  to  the  increased  number  of  calculations. 
There  is  an  optimum  range,  for  each  application,  of  time  step  values  for  the 
integration  scheme. 

1.  Analytical  Solutions 

To  check  the  accuracy  of  the  integration  scheme  as  a  function  of 
time  step,  the  analytical  solution  to  a  second-order  linear  differential 
equation  was  compared  to  a  R-K  4  solution.  The  linear  second-order  actuator 
model  developed  earlier  will  be  used  as  an  example.  Various  time  steps  were 
used  and  the  resultant  average  and  maximum  percent  triors  were  tabulated  and 
plotted.  The  fin  position  and  velocity  (5  and  5)  were  evaluated  using  the 
differential  equations  and  analytical  solutions  presented  earlier.  The  simu¬ 
lation  was  performed  using  a  Zenith  personal  computer  with  an  80286  central 
processor  and  an  80287  numeric  data  processor,  with  all  calculations  done  in 

double  precision.  The  describing  differential  equation  given  earlier  is: 

..  2  o  * 

6  =  ~  _  2£u)q5 

with  a  unit  step  input  (<$c(t)  “  1)  . 

This  equation  will  be  presented  to  the  R-K  4  integrator  to 
arrive  at  a  solution  of  the  first  integral,  velocity  6*  The  analytical  solu¬ 
tion  is: 

5(t)  »  180e  88’4t  sin(n5,2t)  rad/sec 

The  first  integral,  6,  will  again  be  presented  to  the  R-K  4  integrator  to 

calculate  the  second  integral  5,  position  (fin  angle).  The  analytical  solution 

is : 

<5( t )  *  l-1.25e~86’4t  sin(115.2t  +  0.9273)  rad  . 

2.  Simulation  Results 

The  analytical  solutions  and  the  R-K  4  approximations  for  both 
position  and  velocity  were  put  on  the  same  plot  and  are  shown  in  Figures  13 
and  14.  Note  the  visible  errors  in  the  two  plots. 

The  time  step  for  these  plots  was  coarse,  0.01  secs,  and  took  a 
few  seconds  to  complete.  The  percent  errors  at  each  time  interval  were 
averaged  and  are  given  for  position  and  velocity  as  well  as  the  maximum  per¬ 
cent  errors  (see  Table  1). 
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Figure  13.  Plot  of  numerical  and  analytical  solutions  to  fin 
angle  vs.  time,  time  3tep  =  0.01  seconds. 
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TABLE  1.  Table  of  Average  and  Maximum  Percent  Errors, 
h  *  0.01  secs 


Average 

Maximum 

Percent 

Percent 

Position,  6 

0.0511 

3.96 

• 

Velocity,  6 

69.6 

4000 

Next,  a  smaller  time  step  of  0.001  seconds  was  used.  The  errors 
are  given  in  Table  2  and  the  results  are  plotted  in  Figures  15  and  16.  Note 
that  there  are  no  visible  differences  in  the  analytical  and  numerically  inte¬ 
grated  solutions. 


TABLE  2.  Table  of  Average  and  Maximum  Percent  Errors, 
h  =  0.001  secs 


Average 

Maximum 

Percent 

Percent 

Position,  6 

3 . 2xl0~5 

4.7xl0~ 

Velocity,  5 

1.9xl0-3 

0.21 

A  time  step  of  0.0001  seconds  was  used  and  the  plots  are  not 
shown  because  the  errors  were  not  visible,  however  the  errors  are  listed  in 
Table  3. 

TABLE  3.  Table  of  Average  and  Maximum  Percent  Errors, 
h  *  0.0001  secs. 


Average 

Maximum 

Percent 

Percent 

-4 

-3 

Position,  6 

2.4x10 

5.7x10 

• 

_  -7 

-5 

Velocity,  6 

2.7x10 

5.9x10 

Using  a  stop  time  of  0.1  seconds  and  a  resolution  of  0.0005 
seconds,  the  average  percent  errors  for  the  first  and  second  integrals  were 
plotted  as  a  function  of  step  size  in  Figures  17  and  18.  The  maximum  percent 
errors  versus  step  size  were  plotted  in  Figures  19  and  20.  Note  that  ini¬ 
tially,  a  large  step  size  results  in  poor  accuracy.  Smaller  step  sizes  give 
better  accuracy,  as  would  be  expected,  but  too  small  a  step  size  causes  the 
errors  to  increase.  An  optimum  step  size  can  be  chosen  from  these  graphs  or 
from  similar  graphs  from  other  applications. 


The  simulation  was  run  on  a  VAX  11/780  to  investigate  possible 
wordlength  effects  on  the  results;  however,  the  precision  for  the  two  machines 
are  the  same  and  similar  results  were  obtained. 
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Figure  15.  Plot  of  numerical  and  analytical  solutions  to  fin 
angle  vs.  time,  time  step  =0.001  seconds. 
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Figure  16.  Plot  of  numerical  and  analytical  solutions  to  fin 
velocity  vs.  time,  time  step  =  0.001  seconds. 


IV.  CONCLUSION 


The  non-linear  second-order  system  more  closely  modeled  the  actuator’s 
dynamics  and  physical  characteristics  than  did  the  linear  system.  However,  at 
present  we  do  not  have  available  test  data  on  an  actuator's  performance  and 
therefore  cannot  assess  the  validity  of  our  model.  The  Runge-Kutta  integra¬ 
tion  scheme  gave  good  accuracy  and  was  deemed  dependable.  An  optimum  step 
size  of  0.001  second  was  chosen  from  the  percent  error  plots. 
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APPENDIX 

FORTRAN  CODE  FOR 
LINEAR  AND  NONLINEAR  MODELS 


APPENDIX 


The  two  versions  of  the  programs  for  the  fin  actuator  simulation  will  be 
given  here  and  the  subroutines  of  interest  will  be  discussed.  The  linear  and 
non-linear  programs  are  nearly  identical  except  for  the  subroutine  called 
"ACTUAT”,  and  the  analytical  solutions  in  the  main  program  of  the  linear  ver¬ 
sion.  In  both  versions,  two  fins  of  opposite  commanded  deflections  are  simu¬ 
lated,  but  only  the  results  of  the  positive  fin  deflection  are  given  since  the 
negative  fin  deflection  gives  a  mirror  image  response. 

Linear  Model 


In  the  linear  second-order  model,  the  subroutine  "ACTUAT”  contains  the 
differential  equation  that  describes  the  system's  dynamics: 

D£FDD( I)=(0MEGA**2 .DO)*(DEFC(I)-DEF(I)-( 2 .DO*ZETA*DEFD( I) /OMEGA) ) 

This  equation  comes  from  the  differential  equation  developed  earlier. 

The  velocity  state  (DEFD)  is  equated  to  the  variable  X2  for  integration 
purposes . 

The  subroutine  "ACTINT"  initializes  the  actuator  states  and  describes  to 
the  integration  scheme  which  variables  are  integrals  and  integrands. 

Subroutine  "DESOLV"  is  the  Runge-Kutta  integration  scheme.  Variables 
stored  in  the  array  called  "IXDOT"  are  the  variables  to  be  integrated.  The 
corresponding  location  in  array  ”1X"  contains  the  integrated  value.  In  order 
to  perform  integrations  of  orders  higher  than  1,  the  result  of  the  first  in¬ 
tegration  is  stored  as  a  separate  variable  name,  and  this  name  is  submitted  to 
DESOLV  for  integration.  This  continues  until  the  proper  number  of  integra¬ 
tions  have  been  performed.  This  is  the  reason  behind  the  line  of  code: 

DEFD(I)  *  X2(I).  First,  DEFDD  (acceleration)  is  integrated  and  its  integral 
is  called  ”X2".  DEFD  (velocity)  is  assigned  to  the  value  of  X2  and  DEFD  is 
integrated  to  get  DEF  (position). 

Program  "DIGSIM"  is  the  calling  program  and  contains  the  analytical  solu¬ 
tions  and  calculates  the  errors  used  in  plotting.  "C"  arrays  are  used 
throughout  for  communication  between  the  various  subroutines. 

Non-linear  Model 


The  subroutine  "ACTUAT"  in  this  version  has  been  modified  to  include  the 
non-linearities  discussed  earlier.  If  the  absolute  value  of  the  fin  velocity 
is  less  than  1.6  rads/sec,  then  the  rate  feedback  is  zero.  Otherwise  it  is 
( 6  ±  1.6) .  This  is  in  addition  to  the  velocity  state  feedback. 

20 

The  maximum  acceleration  allowed  is  300  rads/sec^.  The  velocity  is 
limited  by  driving  the  acceleration  to  zero  as  the  velocity  approaches  the 
slew  rate  limit  of  5.25  rads/sec.  The  aerodynamic  hinge  moment  is  a  function 
of  the  fin  deflection  and  a  hinge  moment  constant,  ktim. 
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The  position  limiter  sets  the  velocity  and  acceleration  states  to  zero 
and  the  position  state  to  0.436  rads  if  the  fin  is  accelerating  past  the  fin 
stop. 

Note  the  commanded  fin  deflection  is  0.3  rads.  Compared  to  a  unit  step 
(1)  command  input  in  the  linear  model. 

The  rest  of  the  program  for  the  non-linear  model  is  identical  to  the 
linear  program,  with  the  exception  of  the  analytical  solutions,  and  are  not 
given. 
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LINEAR  MODEL 


C 

SUBROUTINE  ACTUAT 
IMPLICIT  REAL*8  <A-H,0-Z) 

COMMON  C (2000) 

DIMENSION  DEF (2)  , DEFD  <2)  ,DEFC(2)  ,X2(2)  , DEFDD (2) 

EQUIVALENCE  (C(  400) , DEF ( 1 ) ) , (C (  401 ) , DEF (2) ) , 

(C<  402 )  , DEFD  < 1 ) )  ,  ( C (  403)  , DEFD (2) )  , 

(C<  404)  , DEFDD <  1) >,  <C(405)  , DEFDD (2) )  , 

<C(  406) ,X2(i> ) , <C<  407) ,X2(2) ) , «C(  1),TIME) 

DATA  OMEGA  /  144. DO  / 

DATA  ZETA  /  .  6D0  / 

DEPC ( 1 )  =  1 . DO 
DEFC ( 2 )  =  —  1 . DO 

C  **  ONLY  TWO  FINS  ARE  SIMULATED  ** 

DO  20  I  =  1,2  « 

C  *SET  DEFD  TO  X2  FOR  INTEGRATION* 

DEFD ( I )  =  X2 ( I ) 

C  *ACTUATOR  DYNAMICS  EQUATION* 

DEFDD  (  I )  =  (0MEGA**2. DO) * (DEFC ( I ) -DEF ( I > -  (2. DO* ZETA* DEFD ( I )  /OMEGA)  ) 

20  CONTINUE 
RETURN 
END 

*-*  THIS  SUBROUTINE  INITIALIZES  THE  ACTUATOR  STATES  ** 

**  AND  SETS  THE  INDEXES  FOR  THE  INTEGRATION  SCHEME  ** 

SUBROUTINE  ACT  I NT 
IMPLICIT  REAL *9  (A-H.O-Z) 

COMMON  C (2000 ) 

COMMON/DEINDX/NDES, IX ( 100) , IXDOTC 100) 

DIMENSION  DEF (2)  , DEFD  <  2 )  , DEFDD (2)  ,X2(2)  ,A<2)  .ADC(2)  ,ADP<2>  , 

ADDC  (2)  ,  ADDF'  (2) 

EQUIVALENCE  (C(  400) ,DEF ( 1 ) ) , (C (  401 ) , DEF ( 2) ) , 

<CC  402 )  ,  DEFD  ( 1 )  )  ,  <  C  <  hO-3)  ,  DEFD  D), 

(C<  404) , DEFDD ( 1 ) )  , (C (405) , DEFDD  <2) >  , 

(C(  406)  ,X2(1) )  ,  <C(  407)  , X2 (2) > 

I  X ( NDES+1 )  =  406 
I X (NDES+2)  =  407 
I  X ( MDE5+3)  =  400 
I  X ( NPES+4 )  =  401 
I XDOT (NDES+1 )  =  404 
I XDOT (NDES+2)  =  405 
I XDO  T (NDES+3)  =  402 
I XDOT 1 NDES+4 )  =  403 
NDES  =  NDES  +  4 
DEF ( 1 )  =  0 . 

DEF (2)  =  0. 

DEFD ( 1 )  =  0. 

DEFD (2)  =  O. 

DEFDD ( 1 )  =  0. 

DEFDD (2)  =  O. 

X2  (  1 )  =  0. 

X2 (2)  =  0. 

RETURN 
END 


A- 3 


n  n  n 


**  MAIN  PROGRAM  ** 

PROGRAM  DIGSIM 
IMPLICIT  REAL*8  (A-H.O-Z) 

COMMON  C ( 2000 ) 

COMMON /DE INDX /NDES , I  X < 100)  , IXDOT (100) 

COMMON/PREQ/  N1 
LOGICAL  LC ( 4000) , QUIT, DATA 
INTEGER  IC (4000) 

DIMENSION  DEF (2) , DEPD (2) 

EQUIVALENCE  (C<i),IC<l>> 

EQUIVALENCE  (C(1),LC(1>) 

EQUIVALENCE  (C(  1>,TJME  ),(C(  24 ) , TSTOP) , ( IC ( 1 20) , NT IME  >, 

(LC (2046) , QU I T ) , (C(19) ,DT> , (C(408> ,Y1) , (C(409) , Y1D) , 
<C<  410) ,DIFFR) , (C<  411) ,DIFFRD) , (C (  400) . DEF ( 1 ) ) . 
(C(  402) , DEFD ( 1 ) ) , (C (  4 12) , DIFAVG) , (C (41 3) , DI FDAVG) 

TSTOP  -  100. 

1  CONTINUE 
QUIT=. FALSE. 

CALL  INPUT (DATA) 

IF (. NOT. DATA)  GO  TO  3 

C - INITIALIZATION 

TIME  =  0. 

NDES  s  o 
NT  I ME  =  0 
Y  l  =0. 

Y 1D=0. 

COUN 1=0 . 

DIFSUM=0. 

DIFDSUM=0. 

CALL  ACTINT 
CALL  ACTUAT 
CALL  OUTPUT 
C - MAIN  LOOP 

2  CONTINUE 
CALL  DESQLV 

C  **  ANALYTICAL  SOLUTION  FOR  COMPARISON  TO  NUMERICAL  ONE  ** 

Y 1 = 1 . ODO-1 . 25D0*DEXP (-86. 4D0*TIME> *DSIN ( 115. 2D0*TIME+. 927295D0) 
Y1D=180. ODO*DEXP (-86. 40D0*T IME) *DSIN( 1 15. 20D0*TIME) 

C  **  DIFFERENCES  IN  THE  ANALYTICAL  AND  NUMERICAL  SOLUTIONS  ** 

DIFFR  = ( DEF ( 1 ) — Y 1 ) /Y 1 
DIFFRD= (DEFD  < 1 ) -Y1D) /Y1D 
COUNT =COUNT  + 1 . 

DIFSUM  =DIFSUM+D IFFR 
D IFDSUM=DIFDSUM+D I FFRD 
DIFAVG  =DIFSUM/COUNT 
D I FD AVG=D I FDSUM/COUNT 
CALL  OUTPUT 
CALL  TERM 

IF ( . NOT.  QUIT)  GO  TO  2 

C - POST  PROCESSING 

CALL  POST  1 
CALL  OUTPTM 
GO  TO  1 

3  CONTINUE 
CALL  F0ST2 
END 
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**  IF  THE  SIMULATION  IS  OVER  ** 

SUBROUTINE  TERM 
IMPLICIT  REAL*8  <A-H,0-Z> 

COMMON  C ( 2000 ) 

EQUIVALENCE  <C(1),LC(1>> 

LOGICAL  LC (4000) , QUIT 

EQUIVALENCE  (C(  1),TIME  ),<C(  24),TST0P  >,<LC< 

IF (TIME. GE. TSTOP) THEN 

PRINT  *, 'SIMULATION  ABORT,  TIME  LIMIT  EXCEEDED: 
QUIT  =  .TRUE. 

END  IF 

RETURN 

END 

BLOCK  DATA  SIM3UBS 
IMPLICIT  REALMS  (A-H,0-Z) 

COMMON  0(2000) 

EQUIVALENCE  <C<  19), DT  ),(C(  170), HDT  ),(C( 

(C(  27) ,TPRSP  >,<C(  40) ,DTPRNT) , (C ( 

(C(  116),  DTPLT  ),.(C(  115),  TPLTO  >,(C( 

(C (  418) , ACTLIM) 

DATA  TSTOP  /  500.  / 

DATA  TPRNTO, TPLTO  /  2*0.0  / 

DATA  TPRSP, TPLTSP  /  2*10.0  / 

DATA  DTFRNT, DTPLT  /  2*  1.0  / 

DATA  DT  /  .00 IDO  / 

DATA  HDT  /  . 0005D0  / 

END 

**  RUNGE  KUTTA  INTEGRATION  SCHEME  ** 

SUBROUTINE  DESOL V 
IMPLICIT  REAL*8  (A-H,0-Z) 

COMMON  C (2000) 

COMMON/DEINDX/NDES, IX (100) , IXDOT ( 100) 

EQUIVALENCE  <C(1),IC(1>) 

REAL*8  K 1 ( 1 00 ) ,  K2 ( 1 00 ) ,  K3 ( 1 00 ) ,  K4 ( 1 00 ) 

INTEGER  IC (4000) 

DIMENSION  XSAVE ( 100) 

EQUIVALENCE  (C(  1),TIME  >,(C(  19), DT  >,(IC(1 

(C(  170) , HDT  ) 

DATA  R6  /. 166666666666666666667D0  / 
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Z046) , QU I T ) 
TIME  =  ' , TIME 

24) , TSTOP  ) . 
.14) , TPRNTO) . 
Ml)  ,  TPLTSP  )  . 


20)  , NT I ME  )  , 


******■**•*■■*■**»****#•****■**#■  »♦*»**»***** 

*  FOURTH  ORDER  RUNGA-KUTTA  METHOD  * 

*.****.****»*-»***■)(• -fr******************** 

C  REVISED  VERSION 

DO  1  It  =  1 , NDES 
J  =  I  X  (  1 1 ) 

XSAVE  (ID  =  C  ( J  ) 

1  CONTINUE 

DO  2  12  =  1 , NDES 
K  =  I XDOT (12) 

K 1 (12)  =  DT*C (K) 

2  CONTINUE 

NTIME  =  NT  I  ME  +  1 
TIME  =  NT IME*HDT 
DO  3  13  =  1 , NDES 
J  =  IX  < 13) 

C(J)  =  XSAVE (13)  +  K1 ( 13) /2. DO 

3  CONTINUE 
CALL  ACTUAT 

DO  4  14  =  1 , NDES 
J  =  IX < 14) 

K  =  I XDOT (14) 

K2 (14)  =  DT  *C ( K ) 

C(J)  =  XSAVE ( 14 )  +  K2(I4)/2.D0 

4  CONTINUE 
CALL  ACTUAT 

DO  5  15  =  1  ,  NDES 
K  =  I XDOT (15) 

K3 (15)  =  DT*C(K> 

5  CONTINUE 

NTIME  =  NTIME  1 
TIME  =  NTIME*HDT 
DO  6  16  =  1 , NDES 
J  =  IX ' 16) 

C(J)  =  XSAVE (16)  +  K3 (16) 

6  CONTINUE 
CALL  ACTUAT 

DO  7  17  =  1 , NDES 
J  =  I  X  (  1 7 ) 

K  =  I XDOT (17) 

K4  (17)  =  DT*C(K> 

C(J)  =  XSAVE ( 17) + (K1 ( 17) +2. DO* <K3( 17) +K3 ( 17) ) +  R4 (  17)  )  /  6.  D< 

7  CONTINUE 
CALL  ACTUAT 
RETURN 

END 
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*******************************************  *******-********-*****  *•*•*■*•*  *■**•**■*•* 

*  OUTPUT  STORES  OUTPUT  DATA  IN  THE  VARIOUS  OUTPUT  FILES  * 

*•#***-*•*  ******■**♦******  *  *****-********■*■********-#*********■*  fr*-***-****-*-*-*** **■■**  * 

SUBROUTINE  OUTPUT 
IMPLICIT  REAL*8  (A-H,0-Z) 

REAL*4  PC (50) 

INTEGER  I C ( 4000 ) 

COMMON  C ( 2000 ) 

EQUIVALENCE  (C(1),IC(1)> 

CHARACTER  NAMES*12 
COMMON / OUTF'SC /  NAMES  (4S)  ' 

COMMON /OUTF'TS/  NS  ,  I NDEXS  ( 48 )  ,  NSPD  ,  INDXSP  (50 ) 

EQUIVALENCE  <C(  1),TIME  >,(C(  26), TPRNT  ),(C(  27),TPRSP  >, 

(C<  40) ,DTPRNT) , (IC(  206), N01  >,(C(  114), TPRNTO ) , 

(C(  116),  DTPLT  ),  <IC<1420  >  ,N02  ),<C(  115),TF'LT0), 

<C<  711)  .TF'LTSP)  ,  (Ct  712)  ,  TPLT  ) 

C  WRITE  PRINT  DATA  TO  LOGICAL  UNIT  61  (SCRATCH  FILE) 

IF(NS.NE.O)  THEN 

IF ( (TIME. BE. TPRNT) . AND. (TIME. LE. TPRSP) )  THEN 
WRITE (61 )  (C ( INDEXS  < I > )  , 1  =  1 ,NS> 

N01  =  NO 1  +  1 

TPRNT  =  TPRNTO  +  FLOAT  (N01  >  *DTF'RNT 
END  I F 
•  END  IF 

C  WRITE  PLOT  DATA  TO  LOGICAL  UNIT  62 
IF (NSPD. NE. 0)  THEN  * 

IF ( (TIME. GE. TPLT) .AND. (TIME. LE. TPLTSP) )  THEN 
DO  1  1=1. NSPD 

1  PC(I)  =  SNGL (C ( INDXSP ( I ) ) ) 

WRITE (62) (PC (I) , 1=1 , NSPD ) 

N02  =  N02  +  1 

TPLT  =  TPLTO  +  FLOAT (N02) *DTFLT 
END  IF 
END  IF 
C 

RETURN 

END 

***•*->(•****•*■*■#■#  *  #***■**# *****■»*■»* **-****■***  *■#■■***# * -4*  *  * •*-*•* *■* *■*■■**♦■*•**■*  ***** 

*  OUTPTM  STORES  THE  MULTI  RUN  DATA  IN  THE  APPROPRIATE  OUTPUT  FILES  * 

*-***  «*♦#»***♦  »*♦*+♦**•»*■»#*■»«*•»«*#■»#■»♦#*#*****  »*»*******♦♦#**■»*****+■»*■*•♦**** 

SUBROUTINE  OUTPTM 
IMPLICIT  REAL*9  (A-H,0-Z> 

REAL*4  PC (50) 

COMMON  C (2000) 

CHARACTER  NAMEM*12 
COMMON /OUTF'MC/  NAMEM  ( 48) 

COMMON /OUTF  M/  NM  ,  I NDEXM  (  48 )  ,  NMPD,  INDXMF'  (50) 

C  WRITE  PRINT  DATA  TO  LOGICAL  UNIT  64  (SCRATCH  FILE) 

I F ( NM . NE . 0 )  THEN 

WRITE (64) (C ( INDEXM ( I ) ) , 1=1 ,NM) 

END  IF 

C  WRITE  MULTI  RUN  PLOT  DATA  TO  LOGICAL  UNIT  65 
IF (NMPD. NE. 0)  THEN 
DO  1  I  =  1 , NMPD 

1  FC  (I)  =  SNGL  (C  ( INDXMF'  (  I  )  )  ) 

WRITE (65)  (PC (  I  )  , 1  =  1 , NMPD) 

END  IF 
C 

RETURN 

END 
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**•#•**•*•*******•*■  **■»■*#*#•»********•***■»**  ************  **********************  ****** 
*  THIS  ROUTINE  WRITES  PRINT  DATA  IN  COLUMN  FORMAT  GENERATED  DURING  A  RUN  * 
***************************  *************  **********  ***■*■*•**■*■*■»****■*  •**  ********* 
SUBROUTINE  POST1 
IMPLICIT  REAL *8  (A-H,0-Z) 

COMMON  C ( 2000 ) 

CHARACTER  NAMES* 12 
COMMON /OUTFSC/  NAMES  (48) 

COMMON /OUTPTS/  NS , INDEXS ( 43 ) , NSPD , I NDXSP (50) 

DIMENSION  X ( 40) 

M=0 

1  CONTINUE 

IF ( (NS. GT. 5*M) . AND. (NS. NE. 0) )  THEN 
M  -  M+l 
I  =  5*  <M— 1 ) +1 
J1  =  5*M 
J  =  J1 

IF(Jl.GT.NS)  J=NS 
REWIND (61 ) 

2  CONTINUE 

WRITE (63, 10) (NAMES (K> ,K=I , J! 

IF ( J 1 . GT. NS)  THEN 

WRITE (63,*)  '  • 

WRITE (63 , *)  '  ' 

END  IF 

DO  3  LINES  =  1,501 

READ (61, END=5 ) (X (K) ,K=1 ,NS) 

WRITE (63 , 20) (X(K) ,K=I,J) 

3  CONTINUE 
da  4  1111  =  lines, 63 

write (63,*)  '  ' 

4  CONTINUE 
GO  TO  2 

5  continue 
do  6  1111  =  lines, 63 

wr  i  te  ( 63  ,  * ) 

6  CONTINUE 
a o  to  1 

END  IF 

CLOSE ( UN I T=6 1 ) 

CLOSE (UNIT=62) 

RETURN 

10  FORMAT (5(2X,A12,2X) //) 

20  FORMAT (5<1X,G14.8> ) 

END 


A-8 


****  I*-*********************************************  ******  *****  ***  **  ********** 

*  THIS  ROUTINE  WRITES  FRINT  DATA  IN  COLUMN  FORMAT  COLLECTED  AFTER  EACH  » 

*  OF  A  MULT-RUN  SET  “■ 

********************************************************************************** 

SUBROUTINE  P0ST2 
IMPLICIT  REAL*8  <A-H,0-Z) 

COMMON  C ( 2000 ) 

CHARACTER  NAMEM* 1 2 
COMMON /OUTFMC/  NAMEM (48) 

COMMON /OUTPM/  NM, INDEXM (48) , NMPD , I NDXMP (50) 

DIMENSION  X ( 48 ) 

M=0 

1  CONTINUE 

I F ( (NM. GT. 5*M) .AND. (NMiNE.O) )  THEN 
M  =  M+l 
I  =  5* (M-l ) +1 
J 1  =  5*M 
J  =  J1 

IF(Jl.GT.NM)  J=NM 
REWIND (64) 

2  CONTINUE 

WRITE (66, 10)  (NAMEM (L)  , L= I , J  > 

IF ( J 1 . GT. NS )  THEN 

WRITE ( 63 , * )  '  ' 

WR I TE ( 63 , * )  '  ' 

END  IF 

DO  3  LINES  =  1,60 

READ ( 64 , END=5 ) ( X ( K ) , K= 1 , NM ) 

WRITE  (66,20)  (  X  ( K )  ,  K.=  I  ,  J  ) 

3  CONTINUE 

dn  4  1111  =  lines, 63 
write (63,*)  '  ' 

4  CONTINUE 
GO  TO  2 

5  continue 

do  6  111!  =  lines, 63 
write (63, *>  '  ' 

6  CONTINUE 
go  to  1 

END  IF 
RETURN 

10  FORMAT (1H1/8(2X.A12,2X) //) 

20  FORMAT (8(1X,G14.8)) 

END 
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SUBROUT  I NE  I NPUT ( DATA ) 

IMPLICIT  REAL*8  <A-H,0-Z> 

COMMON  C  <  2000 ) 

INTEGER  I C (4000) 

LOGICAL  LC (4000) , DATA, LVAL 
D I MENS I  ON  VNAMES ( 50 )  , VNAMEM ( 50 ) 

EQUIVALENCE  (C(1),IC(1>) 

EQUIVALENCE  <C<1),LC(1)> 

EQUIVALENCE  (C (  26 ) , TPRNT ) , ( C (  114),TPRNT0> 

EQUIVALENCE  (IC(206),N01  ) , ( IC ( 1420) ,  N02  ) 

EQUIVALENCE  (C (  712>,TPLT  ) , (C (  115),TPLT0  ) 
character  p  1  t-f  i  1  e  (50)  *8 ,  prnt-f  i  1  (50>*8 
CHARACTER  NAMES* 1 2 , NAMEM * 1 2 , var n am* 1 2 , VNAMS2*8 
CHARACTER  VNAME* 1 4 , VTYPE* 1 , CARD*80 , VNAMES*8 , VNAMEM*8 
COMMON/OUTPSC/  NAMES (48) 

COMMON /OUTPTS/  NS, INDEXS (48) , NSPD , INDXSP (50) 

COMMON/OUTPMC/  NAMEM (48) 

COMMON /OUTF'M/  NM  ,  INDEXM  (  48 )  ,  NMF’D  ,  INDXMP  (50) 

COMMON/FREQ/  N1 
DATA  NFIRST  /  1  / 

DATA  NRUNS/O/ 

data  prnt-fil  /  ’  prntl  *  ,  '  prnt2  ’  ,  '  prnt3  '  ,  '  prnt4  '  ,  '  prnt5  '  ,  '  prn  t6’'  , 

.  ‘ prnt7 ' , 'prnt8' , 'prnt9' , 'prntl O' , 'prntll ' , 'prntl2' , '  prnt 13 ' , 

.  ' prnt 14 ' ,  'prntlS' ,  'prntl6'  ,  '  prnt 1 7  '  ,  'prntl8'  ,  'prntl9'  ,  '  prnt20  '  , 

.  ' prnt21 ' , 'prnt 22' , ' prnt 23 ' , ' prnt 24 ' , ' prnt 25 ' . ' prnt 26 ' ,  crnt27' , 

.  prnt 28' , 'prnt29 ' , 'prnt 30' , 'prnt31 ' , 'prnt 32' , 'prnt 33' , 'prnt34 ' , 
prnt 35' ,  'prnt 36' ,  'prnt 37'  ,  'prnt 38 ' ,  'prnt 39' ,  '  prnt4(.> '  ,  'prnt41  '  , 
.  prnt42' ,  'prnt 43' ,  'prnt 44' ,  ' prnt 45 ' ,  ' prnt 46 ' ,  'prnt47' ,  'prnt 48'  , 
.  ' prnt49 ' , ' prntSO ' / 

data  pit-file  /  '  p  1 1 1  '  ,  '  p  1 12  '  ,  '  pi  t3  '  ,  '  pi  t4  '  ,  '  p  1 15  '  ,  '  pi  t6  '  , 


-Plt7* , 

'pita' , 

pit'*' , 'pltlO' , 'pll 

ill  '  ,  '  pi  1 12  '  ,  '| 

3ltl3 ' , 

'pit  14 ' 

, 'pltlS 

'  ,  ' pi 1 16 

,  'pltl7' , 

'pltlB 

'  ,  '  p  1 1 1 9 

'  ,  '  p 1 120 ' , 

plt21  ' 

, ’ p 1 t22 

', ‘plt23 

,  pi t24 ' , 

'plt25 

'  ,  '  pi t26 

'  ,  ' p 1 127  , 

'  plt2B ' 

, " Pi t29 

'  , ' pi t30 

, ' pi t31  '  , 

'plt32 

'  ,  'pit  33 

'  ,  'Plt34' , 

'  p 1 135 ' 

, ' p 1 t36 

'  ,  ' p 1 137 

, ' pi t38 ' , 

'  pi t39 

'  ,  'pi t40 

'  ,  '  p 1 14 1  '  , 

'  pi t42 ' 

, ' pi t43 

'  ,  'plt44 

, 'pi t45 ' , 

'plt46 

'  ,  ' pi t47 

'  ,  '  pi t48 '  , 

plt49 ' 

, ' pi t50 

'  / 

(n-f  irst 

.  eq.  1 

)  THEN 

NS  =  0 
NM  =  0 
NSPD  =  O 


NMF-D  =  0 
ENDIF 

TPRNT  =  TPRNTO 
TPLT  =  TPLTO 
NO  1  =  0 
NOE  =  0 
N 1  =  0 

NRUNS  =  NRUNS  +  1 
DATA  =  .FALSE. 

nRFN ( 60 , FILE= ' i nput . DAT ' , STATUS= ' OLD ' ' 
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C**** **************** ***********  ****************************  ***  ************* 


c  * 

C  CARD  TYPE  COLUMNS  1-2  * 
C  1- INPUT  DATA  CARD  * 
C  2-PRINT  DATA  CARD  * 
C  3-PLOT  DATA  CARD  * 
C  4-INPUT  DATA  CARD  FOR  LOGICAL  VARIABLES  * 
C  5-DUMP  CARD  * 
C  9-TERMINATOR  CARD  * 
C  10-COMMENT  CARD  * 


*********************************************************************  ****** 


**************************************************************************** 


* 

DECODE  AN  INPUT  DATA  CARD 

» 

* 

DATA  IDENTIFIER  IN  COLUMNS 

4-15 

* 

* 

COMMUNICATION  ARRAY  INDEX  COLUMNS 

16-20 

* 

* 

FORMAT  TYPE  IN  COLUMN 

'T'T 

» 

■* 

FLOATING  POINT  DATA 

26-40 

* 

* 

LOGICAL  AND  INTEGER  DATA 

26-40 

* 

-* 

ECHO  DATA  FLAG  (O-ECHO.  1-NO  ECHO) 

42 

* 

***************************************************************  ************* 
1234  read ( 60 , * , end=999 )  itype 
backspace  60 

i  f  ( itype  .  gt.  10  .or.  itype  .It.  1)  go  to  543 
goto ( 1 00 , 400 , 500 , 543 , 543 , 543 ,543,543, 900 , 1000)  itype 
100  READ  (60,  111)  ITYPE  ,  VNAME  ,  INDX  ,  VTYF'E 


111  FORMAT ' I 2 

,A14, 

14, 

2X , A 1 ) 

if (vtype 

.  eq. 

'F  ' 

.  or  - 

vtype.  eq.  ' -f  '  ) 

go 

to 

101 

i  -f  (vtype 

.  eq. 

'  I  ' 

-  or . 

vtype . eq.  ' i  '  ) 

go 

to 

200 

i f  < vtype 

.  eq. 

'L  ' 

.  or . 

vtype, eq.  ’  1  '  ) 

go 

to 

300 

prints,  'unrecognizable  -format  type' 
go  to  1234 
101  backspace  60 

READ (60, 10)  ITYPE, VNAME, INDX , VTYPE , VALUE , IEDF 
IF ( IEDF  .NE.  1)  THEN 

PRINT*,  ITYPE,  '  VNAME,'  '.INDX,  ',  VTYF'E,'  ', VALUE 
END  IF 

C(INDX)  =  VALUE 

10  FORMAT (I2,A14, I4,2X ,A1 ,E17.8, IX , ID 
go  to  1234  1 

200  backspace  60 

READ ( 60 , 20 )  I TYPE , VNAME , I ND  X , VTYPE , I VALUE , IEDF 
IC(INDX)  =  I VALUE 
IF (IEDF  .NE.  1)  THEM 

PRINT*, ITYPE,  '  '.VNAME,'  ',INDX,'  '.VTYPE,'  ',IC(INDX) 
END  IF 

20  FORMAT  < I2,A14, I4,2X,A1 , 10X , 1 7, 1 X , 1 1 ) 

GO  TO  1234 
300  backspace  60 

READ ( 60 , 30 )  I TYPE , VNAME , I ND  X , VTYPE , LVAL , IEDF 

LC(INDX)  =  LVAL 

IF ( IEDF  .NE.  1)  THEN 

PRINT*.  ITYPE,  '  '.VNAME,'  '.INDX,'  '.VTYF'E,'  '  ,  LC  ( INDX  ) 
END  IF 

30  FORMAT ( 12 , A1 4 , I  4 , 2X , A1 , 2X , L7 , 9X , 1 1 ) 

GO  TO  1234 
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C***********************************************************************"1'-*'*  * 
C  DECODE  A  PRINT  DATA  CARD  * 

C  PRINT  HEADER  IN  COLUMNS  4-15  * 

C  COMMUNICATIONS  ARRAY  INDEX  IN  COLUMNS  16-20  * 

C  ,PR I  NT  DATA  COLLECTION  FLAG  IN  COLUMN  23  * 

C  MAXIMUM  OF  48  PRINT  DATA  CARDS  * 

C************ ********************** ***************************  ************** 
400  PZmD  < 60 , 40 )  IT YFE  ,  var nc ..n  ,  I  I NDX  ,  MF  FLO 
Backspace  60 

i-f(mpflg  .  eq.  1)  go  to  401 

NS  =  NS  M 

IF (NS  .GT.  48)  THEN 

PRINT*. 'TOO  MANY  PRINT  VARIABLES  ' 
go  to  1234 
endi  f 

READ  (60, 40)  I  TYPE  .NAMES  (NS)  ,INDEXS(NS)  ,  MF'FLG 
40  FORMAT (12, 1X,A12,1X,I4,2X,I1) 

PRINT*,  ITYPE,  '  NAMES  (NS)  ,  '  '  ,  INDEXS  (NS)  ,  '  '  ,  MF'FLG 

GO  TO  1234 
401  NM  =  NM  +  1 

IF (NM  .GT.  48)  THEN 

PRINT*,  TOO  MANY  PRINT  VARIABLES  ' 
go  to  1234 
ENDIF 

READ (60,40)  ITYPE , NAMEM (NM) , INDEXM (NM) , MPFLG 
PRINT*, ITYPE,  '  ', NAMEM (NM)  ,  '  ', INDEXM (NM)  ,  '  ', MPFLG 

GO  TO  1234 

C**********  *******************************************  *****  ***************** 


C  DECODE  PLOT  DATA  CARDS  * 
C  PLOT  LABEL  -4-15  * 
C  C  ARRAY  INDEX  COLUMNS  16-20  * 
C  PLOT  DATA  COLLECTION  FLAG  23  * 
C  THE  PLOT  DATA  COLLECTION  FLAG  INDICATES  WHETHER  PLOT  DATA  13  * 
C  TO  BE  COLLECTED  THROUGHOUT  EXECUTION  OF  A  SIMULATION  RUN  OR  * 
C  IS  TO  BE  COLLECTED  ONLY  AFTER  EXECUTION  HAS  BEEN  COMPLETED  * 
C  MAXIMUM  OF  50  PLOT  DATA  CARDS  * 


C **************** ************************  *****************************  *  ***** 

500  READ (60,60)  ITYPE , VNAMS2 , I INDXSP , MPFLG 
backspace  60 

if(mpflg  .eq.  1)  go  to  501 
NSF'D  =  NSF'D  +  1 
IF (NSPD  . GT.  50)  THEN 

PRINT*, ’TOO  MANY  PLOT  VARIABLES  ' 
go  to  1234 
ENDIF 

READ (60, 60)  I  TYPE , VNAMES ( NSPD )  , INDXSP (NSPD)  , MPFLG 
60  FORMAT (12, A8 ,6X,I4,2X,I1) 

PRINT*,  ITYPE,  '  VNAMES  (NSF'D)  ,  '  ',  INDXSP  (NSPD)  ,  '  ',  MPFLG 

GO  TO  1234 

501  NMF’D  =  NMPD  +  1 

IF  (NMF'D  .GT.  50)  THEN 

FRINT * ,  ' TOO  MANY  PLOT  VARIABLES 
go  to  1234 
ENDIF 

READ  (60,60)  I  TYPE  ,  VNAMEM  ( NMPD )  ,  INDXMP  (NMPD)  .MF'FLG 
PRINT*,  ITYPE,  *  ',  VNAMEM  (NMPD)  ,  '  INDXMP  (NMPD)  ,  '  ', MF'FLG 

GO  TO  1234 
543  read (60,90) card 

print*, 'card  type  not  found' 
print*, card 
print*, ' 
go  to  1234 

1000  read (60, 90) card 
print*, card 
go  to  1234 
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c  *  #■»■»****■*  *  *  *  »*****»*»**4****#4  *■*•*•*■■****#**#**  **■-**  »***♦**♦»»  «"**«*-**»  ********* 

C  TERMINATOR  CARD  MUST  FOLLOW  DATA  DECT  FOR  EACH  RUN  * 

C  RUN  TITLE  IS  CONTAINED  IN  COLUMNS  4-23  * 

C******************************************************************  ********* 
900  READ (60,90)  CARD 
PR I NT*, CARD 
90  FORMAT (A) 

IF (NS  . NE.  0)  THEN 

OPEN (61, FORM= ' UNFORMATTED ' , ST ATUS= ' SCRATCH ' ) 

OPEN ( 63 , F ILE=prnt f i 1  (NRUNS)  , ST  ATUS= ' UNKNOWN ' ) 

END  I  F 

I F ( ( NM  .NE.  0)  .AND.  (NFIRST  . EQ.  1>)  THEN 
OPEN (64 ,FORM= 'UNFORMATTED ' , STATUS= ' SCRATCH ' ) 

OPEN (66 , FILE= ' FRINTM. DAT ' , STATUS= ' UNKNOWN ' ) 

END  I F 

IF (NSRD  .NE.  0)  THEN 

OPEN ( 62 , F ILE=p 1 1  f i 1 e ( nruns)  , FORM= ' UNFORMATTED  ) 

WRI TE (62)  NSPD.NSFD 

WRITE (62)  (VNAMES (III)  , I  I  1  =  1 , NSPD) 

END  IF 

IF<  (NFIRST  .EG).  1)  .AND.  (NMPD  .  NE .  0)  )  THEN 
OPEN ( 65 , F I LE= ' PLTM . DAT ' , FQRM= ' UNFORMATTED ' ) 


WRITE (65)  NMF 
WRITE (65)  (VN 
END  IF 
NFIRST  =  0 
DATA  =  .TRUE. 

999  CONTINUE 
RETURN 
END 


NMPD, NMPD 
(VNAMEM (III) 


I  I  1  =  1 , NMPD ) 
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NONLINEAR  MODEL 


SUBROUTINE  ACTUAT 
IMPLICIT  REAL*S  (A-H,0-Z) 

COMMON  C (2000 ) 

DIMENSION  DEF (2)  , DEFD (2)  ,DEFC<2)  , X2(2)  .ACCRQ(2)  , RATEF3 ( 2)  , 
ACCEL IM1 (2) , ACCEL I M2 (2) ,HM(2> , DEFDD (2) 

EOu  I  v  Ai_2NCE  <  C  (  400  )  ,  DEr  1 1 )  >  ,  <  C  4u  1  >  ,  DEF  ( 2  >  /  , 

(C(  402) , DEFD (1) ), (CC  403) , DEFD (2) ) , 

(C(  404) , DEFDD ( 1 ) ) , (C ' 405) . DEFDD (2) ) , 

(C(  406)  ,X2(1) >  ,  <C(  407)  ,X2(2> )  ,  CC  <  1),TIME> 
DATA  RATELM  /  5.25D0  / 

DA’ A  OMEGA  /  144. DO  / 

DATA  ZETA  /  . 6D0  / 

DATA  32  /  300. DO  / 

DATA  31  ,■  TS7.3D0  / 

DATA  K:HM  /  0 .  DO  / 

DEFC ( 1 )  =  . 30 DO 
EEFC (2)  =  30D0 

***  ONLY  TWO  ACTUATORS  ARE  SIMULATED  *** 

DG  20  I  =  1.2 

♦SET  DEFD  TO  X2  FOR  INTEGRATION* 

DEFD  (I )  =  X2  (  I ) 

♦CALCULATE  DEADBAND  IN  RATE  FEEDBACK* 

RATEFB(I)  =  0 . DO 

IF  (DEFD (I ) . 3T. 1 . oDO)  RATEFB ( I )  =  (DEFD ( I ) -1 . 6D0) /20. DO 
IF  (DEFD (I) . LT.-1.6D0)  RATEFB ( I )  =  (DEFD ( I ) +1 . 6D0) /20. DO 
♦ACTUATOR  DYNAMICS  EQUATION* 

ACCRQ ( I )  =  ( CMEGA**2 . DO ) * ( DEFC ( I ) -DEF  < I ) 

- (2. DO*ZETA*DEFD ( I ) /OMEGA) -RATEFB ( I ) ) 

♦CALCULATE  SPEED-TORQUE  LIMITS* 

ACCEL  I Ml ( I )  =  31* < 1 . DO-DEFD ( I ) /RATELM) 

IF  (ACCELIM1  (I)  .GT.G2)  ACCELIMKI)  =  G2 
ACCEL I M2 ( I )  =-Gl* ( 1 . DO+DEFD ( I ) /RATELM) 

IF  ( ACCE- I M2 ( I ) . LT . -G2 )  ACCEL I M2 ( I )  =  -G2 
IF  (ACCRQC I)  .3T.ACCELIM1  •  I)  )  ACCRQ  ( I )  =  ACCELIMKI) 

IF  (ACCRQC)  .LT.  ACCELIM2(I)  )  ACCRQC)  =ACCELIM2(I) 
♦CALCULATE  HINGE  MOMENTS* 

HM (I )  =  DEF ( I > *KHM 

♦CHECK  FOR  POSITION  LIMIT* 

IF  (DEF (I ) . 3E. 0. 436D0  .AND.  DEFDD (I > .ST.O.DO)  THEN 
DEFDD ( I )  =  0. DO 
DEFD  C)  =  0.  DO 

DEF  ( I )  —  (.).  436D0 
END  IF 

IF  (DEFC)  .LE.-0.436D0  .AND.  DEFDD  (  I )  .  LT.  0.  DO)  THEN 
DEFDD (I)  =  0 . DO 
DEFD (I)  =  0. DO 
DEF  .1)  —  — . 436D0 
END  I F 

DEFD  ( I )  =  X2C) 

DEFDD (I)  =  ACCRQ C)-HM(I) 

CONTINUE 

RETURN 

END 


**  THIS  SUBROUTINE  INITIALIZES  THE  ACTUATOR  STATES  ** 

**  AND  SETS  THE  INDEXES  FOR  THE  INTEGRATION  SCHEME  ** 

-UBRCUTINE  ACT I. NT 
IMPLICIT  REAi_*8  (A-H.C-Z) 

COMMON  0(2000) 

CGM-CN/DEINDX/NDE3. IX ( 100) . I XDCT ( 1 00 / 

DIMENSION  DEF  2)  ,  DEFD  (  2  )  .DEFDD  (2)  .X2’2>  ,A*.2>  ,  ADC  (2)  .  ADR  '  2)  . 
ADDC  2' ,  ADDR (2 ■ 
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